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1. Introduction 



The basic theory of Fresnel diffraction at plane aper- 
tures was developed long ago [1,2] and is summarized in 
textbooks [3-6]. For apertures bounded by straight lines 
(rectangle, slit, half plane), the standard textbook solu- 
tion in terms of complex Fresnel integrals' is based on 
special but poorly documented transformations of coor- 
dinates. It is shown in this paper that such transforma- 
tions cannot be performed accurately for apertures irra- 
diated by arbitrarily located point sources. In the past, 
the numerical evaluation of the complex Fresnel inte- 
grals themselves has also been a problem, and thus pre- 
vious discussions were confined to simplified special 
cases. Computational details were omitted and semi- 
quantitative methods (Cornu spiral) were used to de- 
scribe the nature of diffraction at rectangular apertures. 

In the case of circular apertures, the rigorous solution 
involves Lommel functions of two variables, which are 



^ In this paper it is necessary to distinguish between the "Fresnel 
diffraction integral" Uy{P) defined by Eqs. (3a-c) of this paper and 
the "complex Fresnel integral" F{s) defined by Eq. (8). 



defined as series expansions in Bessel functions and 
previously had to be evaluated by tedious manual calcu- 
lations or approximations. For the most part, these ap- 
proaches have been rendered obsolete by modern com- 
puter software. However, approximative methods are 
still useful for work on personal computers which in- 
volves large values of the configuration parameter u 
defined by Eq. (17a) of this paper. It is shown here that 
a previously used approximation by Focke [7] is inade- 
quate for this purpose on account of its poor accuracy, 
but that an older approximation by Schwarzschild [8] 
gives excellent results. 

As algorithms for the computation of Fresnel diffrac- 
tion patterns on a personal computer have not been 
published, a compilation of such algorithms is presented 
in this paper. The underlying theory is stated for off-axis 
source points, so that the results can be applied to ex- 
tended sources. For both types of aperture, the closed 
solutions obtained are paraxial approximations. 
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2. The Fresnel Diffraction Integral 

The scalar wave function U(P) associated with dif- 
fraction at a plane aperture is customarily expressed by 
one of the Rayleigh-Sommerfeld integrals^, 



UUP)=-^jdQ^' 



KPoQ+QP) 



2lT 



., 1 \ dPoQ ,, , 



t/UP)=-^/dGg9(i.-J^)^, (lb) 



or, alternatively, by the Kirchhoff integral, 

UK(P) = \(Uh+Uis). 



(Ic) 



Here, as indicated in Figs. 1, 2, and 4, Pq is the location 
of a point source emitting a monochromatic spherical 
wave of amplitude Agph, circular wave number k = lir/X , 
2 is a point in the aperture, dQ is the surface element at 
2, w is the aperture normal pointing away from the 
source, and P is the point of observation. 

In the Fresnel approximation the points Pq and P are 
located at finite distances which are large compared to 
the wavelength of light and the dimensions of the aper- 
ture. Therefore, it is assumed that 



1 ,1 , 



(2a) 



and that the distances PoQ and QP and their normal 
derivatives do not vary appreciably inside the aperture. 
Thus they can be replaced, except in the rapidly oscillat- 
ing exponential function in the integrand, by their values 
at an arbitrarily chosen reference point O inside the 
aperture. Under these conditions, the first Rayleigh- 
Sommerfeld integral Eq. (la) may be written in the form 



UUP) - 



_ iMsph dPpO e^ 



ik(PoO+OP) 



2tt dn PoO'OP 



\ 



dQe 



ifcA(e) 



where 



KQ) = {PoQ + QP)-{PoO + OP\ 



(2b) 



(2c) 



is a small quantity that can be expressed in approximate 
form. It is well known that, with 



^ A subscripted notation for wave amplitudes is used to avoid confu- 
sion between quantities which differ in physical significance and di- 
mension. For example, the squared amplitude lAsphl^ of a spherical 
wave denotes a radiant intensity whereas the squared amplitude 
lApiancI^ of a plane wave denotes an irradiance. 



O = (0, 0, 0), Po = U, yo, Zo), Q = it V, 0), P = {x,y,z) 

(2d) 

expressed in cartesian coordinates, the required approx- 
imation for A(G) is 



A(G) =-[(/- /o)f + (m - mo)7,] + ^ He + V') 



(2e) 

where 6(^,17) is the residual error when terms of third 
and higher order in ^ and j] are neglected. Here, 



/o= ,wio= ^ l = - ^ m = - (2t) 

To Vq r r 



are the first and second direction cosines of the vectors 
Poland OP, and 



To = Vxo + jo + zo, r = Vx^+j^ + z^ 



(2g) 



are the distances PqO and OP . The corresponding value 
of the normal derivative^ in Eq. (2b) is 



-T^ — = ~ ^^= =~ COS^o- 



(2h) 



We now have 



UUP) ■ 



iMsphCos^o 
li^rov 



JA(ro+r) I ^ 



dGe^ 



ifcA(e) 



(2i) 



The corresponding forms of UUP) ^nd U^{P) are es- 
sentially the same, except that — cos^o is replaced by 
cos^and V2(— cos^o + cos^), respectively, where ^is the 
colatitude of the point of observation P. Because the 
diffracted light is confined to a narrow angular range 
about the central direction PqQ unless the aperture di- 
mensions are extraordinarily small, these differences 
may be judged insignificant. As the Rayleigh-Sommer- 
feld solutions pertain to the respective cases of/?- and 
5 -polarization of the incident light, this implies that 
Fresnel diffraction is independent of polarization. The 
Kirchhoff solution has no definable meaning as far as 
polarization is concerned, but turns out to be equivalent 
to the Rayleigh-Sommerfeld solutions in the Fresnel 



^ The angle 6*0 should not be confused with the angle tt — 6*0 indicated 
in Figs. 1 and 4. 6 and 6*0 are colatitudes which are measured clockwise 
from the positive z -axis. In this paper, 6*0 is assumed to be on the order 
of 77, so that cos 6^ is on the order of —1. 
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approximation. A further solution, the Maggi-Rubinow- 
icz transformation of Kirchhoff's integral [3,4], is not 
suitable for computations of Fresnel diffraction patterns 
because it is singular at the boundary of the geometrical 
shadow. 

In this paper, Eq. (2i) will be regarded as the basic 
form of the Fresnel diffraction integral and will be writ- 
ten as 



(l^-\-m7]y are much smaller than (^^ + 17^), and then 
one finds'^ 



e(f,T?) = 






ie+vy 



(4b) 



Hence, the magnitude of 6(^,17) relative to the quadratic 
term of Eq. (2e) may be estimated as 



Up(P)=- Uo(P)cosOoa^(PX 



where 



Uo(P)=A 



sph 



Jk(ro+r) 



ro + r 



: VEo(Py 



ik(ro+r) 



(3a) 



(3b) 



is the geometrical field at the point of observation P 
according to Huygens' principle. 



ap(P) 



^ ik(ro-\-r) f 
lirror J 



dQd 



}kMQ) 



(3c) 



is the modification of the geometrical field by diffrac- 
tion, Eq(P) is the normally incident geometrical irradi- 
ance at P, and — cos^o is the inclination factor according 
to Lambert's law. 

The third- and fourth-order terms neglected in Eq. 
(2e) are 



€(^,V)=- ^AUo^^movMe^v") - Uo^^movf]} 



+ 2^2 m^mr^Me + V') - d^^mr^f]} 



- 8^3 ae + v'f - (io^^movf[6(e + v') 



2roH6(g,7?)l (rj^r'Xe^v') . ^. 



2 

max 



(ro+rXe+rj^) 4rSr\ro+r) ^ {rf ' 



(4c) 



where ^^ax is the maximum value of V^^ + 7]^ (e.g., the 
radius of a circular aperture) and (r) is an average of ro 
and r. Accordingly, the relative error in A(2) is in- 
versely proportional to the square of the relative distance 
{r)/qmax- At a distance of ten aperture dimensions, it is 
on the order of 1 %. 



3. Rectangular Aperture 
3.1 General Theory (Fig. 1) 

When applying the above equations to a rectangular 
aperture of width 2w and height 2h , it is customary to 
transform the global cartesian coordinates (x, j, z) as- 
sumed in Sec. 2 into local coordinates (x', j', z') which 
depend on the locations of the points Pq and P and are 
chosen so that Eq. (3c) is separated into a product of 
independent Fresnel integrals in ^ and 7]. That is. 



a^(P) ' 



jv^-^j 



d7]d 



ibrf^ 



(5) 



The first step in this transformation is to place the origin 
of the local coordinates at the point M where the straight 
line PqP intersects the aperture plane: 



5(/o^+moT7)']} 



^{(e^Vy-(l^^rn7^f[6(e^v') 



M = (xm, Jm, 0), Xm = , Jm = - — _ . (6a) 

Z Zo Z Zo 



This gives 



-Sd^^mr^f]}. 



(4a) 



This shows that the error introduced by neglecting this 
term depends in a complicated manner on the geometri- 
cal parameters involved. Accordingly, it is difficult to 
assess its magnitude without considering specific cases. 
However, a few general comments are in order. Under 
ordinary circumstances, the direction cosines Zq, /, ... are 
small compared with unity, so that (/o^+^oi?)^ and 



, _ Xq Xm _ jr _ ^ ^M 

Lq — ; — I — ; 

ro r 



ro r 



ro = 



Zo 



cosOm 



,r = 



cosOm ' 



(6b) 



"^ The direction cosines are identically equal to zero, and Eq. (4b) holds 
exactly, for the rectangular aperture discussed in Sec. 3.1.2. See Eqs. 
(9a-d). 



499 



Volume 103, Number 5, September-October 1998 

Journal of Research of the National Institute of Standards and Technology 



Om being the angle indicated in Fig. 1 . The linear term 
of Eq. (2e) now vanishes and we have 



HQ) = ^A(^-^Mf^(v-yMf 



Uo'i^ - Xm) + mo'(7] - yM)f 



P =7T 



-ZZo 



ro' + r' (z - Zo)cos^ 



(6c) 
(6d) 



This result can be used in two ways to derive a final 
result. 



3.1.1 Paraxial Approximation 

As mentioned in deriving Eq. (4b) above, the direc- 
tion cosines /o' and mo' will be small if the points Po and 
P are close to the z-axis of Fig. 1. To a first-order 
approximation in Om we have I'o^, m'^ « 1, so that the 
third term of Eq. (6c) can be omitted and Eq. (3c) leads 
directly to 



. , pw rh 



ik(r)-yMr/2p' 



= -- [F(sd - F(s-)][F(h) - F(t-)l (7a) 



where 



k k 
s± = \ ;(± w - Xm), t± = a/ ;(± h - Jm), (7b) 

\ 7Tp \ 7Tp 



and 



F(s) = C(s)-^iS(s) 



= j do-i 



jra^/2 



(8) 



is the complex Fresnel integral. 

3.1.2 Coordinate Transformation for Off-Axis 
Sources 

According to textbooks, a rotation of coordinates may 
be necessary when the direction cosines /q' and mo' 

I 



i 

2w 


1 


I 


2h/ 9 

/ o 


\ 


/ 


// 


' m/ q/ ' " 


x-9„ 


\ 1 
■ 1 





Fig. 1. Notation for rectangular apertures. 

are too large to justify the paraxial approximation of 
Sec. 3.1.1, a rotation of coordinates may be necessary. 
The usual recommendation [3, 5] is to place the new 
x'-axis along the projection of the line PqP onto the 
aperture plane. This gives mo' = 0, so that a^^iP) does 
indeed assume the form stipulated by Eq. (5). However, 
the x'- and j'-axes so defined are not parallel to the 
edges of the aperture, and consequently the two inte- 
grals are not separable because the limits of the ^-inte- 
gral depend on 17, and vice versa. 

To overcome this difficulty, a different transformation 
is attempted in the following: The z'-axis is placed in the 
direction of the unit vector along the line P^P , so that 
/o' = niQ = and Eq. (5) is again satisfied. The x'-axis is 
chosen so that its projection onto the aperture plane is 
parallel to the ^-direction. The j'-axis is defined in the 
usual manner as y' = z' X x'. Thus, 



fx\ / l/W 

y' = - sin</>MCOs</>Msin^Mtan^M/W WcosOm 
\z7 \ cos</>Msin^M sin</>AfSin^M 



- cos</>Mtan^M/W\ /x - Xm\ 

- sin(/>Msin^M/W 11 j - Jm 

COS^M /\ z j 



(9a) 
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where cfy^ and Om are the longitudes and colatitudes of Pq 
and P with respect to M, 



^ ^0 Z Zo 



and 



W= Vl + cos^f^A/tan^^M. 
Accordingly, for z = 0, 



^'- w (^^ ^^)' 



(9c) 



r]'= VKcosft, 



sin(^MCQS<^Mtan^^M ^^ x 



iv^ 



+ (l? - Jm) 



(9d) 



Equation (9d) shows that 7]' is still not independent of ^. 
This was to be expected as it is not possible to rotate the 
z-axis and have orthogonal x- and j -axes which are both 
aligned with the aperture edges. Accordingly, the sepa- 
ration of integration limits is not complete unless the 
first term in the above expression for 7]' is omitted — a 
first-order approximation in Om- Then, 



^ = — (^-Xm), r]'=WcosOM(r] 



d^'dT7'= cos^Md^dr?. 



■ Jm), 



(10a) 



and 



a,(P) = 



ikCOS Om ^^^i^(^-jCM)2/2W2p'2 I ^ ifcw2cos2^r,-3^M)^/2p'2 



27TP' 



j-h 



= - - [F(sJW) - F(sJW)][F(WcosOMtd 



- F(WcosOMt-)l 



(10b) 



where s+ and s.. are the same as in Eq. (7b), above. It 
should be noted that this result differs from the paraxial 
approximation only by the factors l/W and VKcos^m in 
the arguments of the Fresnel integrals. Within the above 
approximation for 17' these factors are equal to unity, and 
thus Eq. (10b) appears to be no improvement over the 
paraxial approximation Eq. (7b) of Sec. 3.1.1. It follows 
that, for rectangular apertures, the coordinate transfor- 
mations recommended in Refs. [3] and [5] are superflu- 



ous. To higher than first order in Om, the ^- and 17 -inte- 
grals remain inseparable and a closed solution for ap(P) 
is not possible. 

3.1.3 Evaluation of Fresnel Cosine and Sine 
Integrals 

The use of Eqs. (7a) and (10b) is straightforward. An 
example is given in Sec. 3.2, below. It should be remem- 
bered that the variation of a^ (P) with P is implicit in 
Eqs. (7b) and (10b), in that x^, Jm* 4^m and Om depend on 
the location of P. It should also be borne in mind that the 
point M of Fig. 1 will be outside the aperture when P 
lies in the geometric shadow. This can lead to values of 
^ and 7] larger than assumed in Eq. (3a). For this reason, 
the computation of ap (P) based on Eq. (7b) or (10b) 
must not be carried too far into the shadow region. 

The only problem that may be encountered on a per- 
sonal computer is that the Fresnel cosine and sine inte- 
grals defined by Eq. (8) are not usually included in 
standard software packages. For modest accuracy re- 
quirements, they can be computed from the equations 
quoted in Ref. [9], 

C(5) = ^+/(5)sin(^j - g(s)cos\^l 

C(- s) = -C(sX (11a) 

S(s) = ^-f(s)cos\^] - g(s)sm[^l 

S(- s)=-S(sX (lib) 

1+0.926 5 



f(s) = 



2+1.792^ + 3.104 5' 



+ 6(5), 5^0, (lie) 



^^^^ = 2 + 4.142. + 3.492.^6.67.^^^^^^'^^^' 

(lid) 



where le(.)l ^ 2 X 10 ^. Accordingly, the following 
simple algorithm may be used: 

1. Define .. 

2. Let .'=1.1. 

3. Calculate /(.'), g(s') from Eq. (llc,d). 

4. Calculate C(.'), S(s') from Eq. (lla,b). 

5. If . < let C(.) = - C(.'), S(s) = - S(s'). Else, let 
C(s) = C(s'XS(s) = S(s'). 

If better accuracy is desired, this algorithm can be im- 
proved by using the method described in Ref. [10]. 
Alternatively, software for computing C(.) and 5'(.) in 
Fortran or C can be downloaded [11, 12]. 



501 



Volume 103, Number 5, September-October 1998 

Journal of Research of the National Institute of Standards and Technology 



3.2 Application to Slits (Fig. 2) 

The rectangular aperture discussed so far is trans- 
formed into a slit of width 2w on setting h = ^ in Eq. 
(7b). ^ It may then also be assumed that the source is a 
long luminous line which is parallel to the slit and passes 
through the point Pq in Fig. 1, so that it will suffice to 
compute the diffraction pattern in the xz -plane shown in 
Fig. 2. With these assumptions we have ^+ = ±00^ so that 
F(U) = ± V2 (1 -F i), [F(h) - F(t-)] = 1 + i, and Eq. 
(10b) is reduced to 



^M — 0, Um — Oc 



a^(P) = 



1 



■ [F(sd - F(s-)] 



1 -i 



{[C(s,) - C(s-)] + i[S(s,) - S(s-)]}, (12a) 



with s as defined by Eq. (7b) but assuming 
yQ = y = y^ = Q SO that Eq. (9b) is simplified to 



6m = arctan , x = Xm + ztanft^. (12b) 

Zo 



As the diffraction pattern is centered at and symmetrical 
about the geometrical source image C shown in Fig. 2, 
where 




: arctan — , x = 

Zo 



■Xc' 



_ XqZ 

Zo 



(12c) 



it will also suffice to compute it for positive value values 
of Xm, only. The computation is typically carried to a 
maximum value of 9m beyond the shadow boundary S , 
the latter being given by 



Xm = w, 0m= Os = arctan — ^ x = xs = w -\- ztanft. 



Zo 



(12d) 



Accordingly, the following procedure may be used to 
evaluate the dependence of a^iP) on x: 

1. Define a maximum (xM)max and a step size Ax^ for 
Xm- 

2. Let Xm = 0. 

3. Compute 6m and x from Eq. (12b) and s± from Eq. 
(7b). Use the algorithm of Sec. 3.1 to find C(s±) and 
S(s±). Compute ap(P) from Eq. (12a). 

4. LetxM = Xm + Axm. If Xm < (xM)max, go to Stcp 3. Elsc, 
stop. 

A typical diffraction pattern computed in this manner is 
shown in Fig. 3. The numerical parameters chosen for 
this particular example are listed in the figure caption 
and were taken from an experiment described by 
Fresnel.^ 



6t 




Fig. 3. Relative irradiance, la^iu, v)P vs v/«, for a slit, w = 1 mm, 
Zq=- 2.507 m, z = 1.140 m, A = 639 nm. 



Fig. 2. Notation for slits. 



^ This assumption seems to violate the condition that the values of rj 
in Eq. (2e) must be small, but is justifiable on account of Fresnel's 
zone construction. Because the field at P is not affected by zones 
located at large distances from the line PqP in Fig. 1 , the 17 -integration 
can be extended to infinity without introducing an error. 



^Reference [1], pp. 117 and 128. At first glance, the variation of 
\aFiP)\^ shown in Fig. 3 does not appear to match the results described 
by Fresnel. However, the agreement is satisfactory when the data are 
replotted on a logarithmic scale to simulate Fresnel's visual readings. 



502 



Volume 103, Number 5, September-October 1998 

Journal of Research of the National Institute of Standards and Technology 



4. Circular Aperture (Fig. 4) 
4.1 General Theory 

In evaluating the Fresnel diffraction integral Eqs. 
(3a-c) for a circular aperture with diameter 2a it is 
convenient to use spherical coordinates centered at the 
aperture center O, so that 

Po = (xo, Jo, Zo) = ro(cos(^o sin^o, sin<^ sin^o, cos^o) 

= - ro(lo, mo, TioX (13a) 

Q = (^, V, 0) = ^(cos;^, sin;^, 0), (13b) 

P = (x, J, z) = r(cos(^ sin^, s'mcf) sinO, cos6) 

= r(l, m, n). (13c) 

In these coordinates, the path difference A(2) defined 
by Eq. (2e) can be evaluated as follows. 



2a 




Let C = r(lo, ruo, Hq) be the geometrical image of Pq at 
the distance r from the aperture center, so that the posi- 
tion of P relative to C will be given by the vector 
CP = r(l — lo, m — mo, n — rio). Let F be the foot of the 
perpendicular from C onto the xy -plane, let c = FP , and 
define 

FP = r(l — /q, m — mo, 0) = c(cosy, siny, 0), (14a) 

where 

c = rVil - kf + (m - mo)' 

= rVsin'^+ sin'^o + 2sin^sin^ocos((^ - (^), (14b) 
m — mo sin(^sin^+ sint^osin^o 



tany = 



I — lo cos(^sin^+ cos(^osin^o 



(14c) 



Accordingly, the linear term of Eq. (2e) can be ex- 
pressed in the form 



(/ - /o)^+ (m - mo)7] = ^FP-OQ 



qc qc 

■ -^- (cos;^cosy + sin;^siny) = ■^-~ cos(;^ — y). (14d) 



In the quadratic term of Eq. (2e), we have 

= q\l — sin'^o(cos(^ocos;^+ sint^osin;^)'] 

= q\l - sin'^ocos'(;^f - c^o)] (15a) 

and, likewise, 

ie + V') - il^+rni^f = q\l " sin'^cos'(;^ - <!>)]. 

(15b) 

In the following, it will be assumed that the points Po 
and P are close to the z-axis so that sin'^o and sin'^ are 
negligibly small compared to unity. In this paraxial ap- 
proximation one obtains 

^^[(e^v')-Go^^movf] 



+ ^[(^'+77')-(/^+m77)'] 



2l T0 + ^ 

2ror 



(15c) 



On substitution of Eqs. (14d) and (15c) into Eq. (2e) we 
have 



Fig. 4. Notation for circular apertures. 
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A(G): 






(16a) 



and hence Eq. (3c) is reduced to 

.„. ik(ro-\-r) T, ifr^s±^«2 f^, _.u9£ 



cos(^-7) 



= _ JM^ ^\qqj,(k^y^^\ (16b) 



where the integral over x was evaluated as lirJoikqc/r) 
[9]. On substituting 



a ka (ro-\- r) kac 

p = -, u = , v = 

a r^r r 



(17a) 



4,2 Lommel's Solution 

Lommel [2] evaluated the integral Eq. (17b) in the 
form' 






dpp/o(vp)ei«^ = i VUu, v) + iM(i/, v)], (18a) 



so that 



q;f(m, v) = aiiu, v) = - M(m, v) - \^l^{u, v). (18b) 



In this notation, 



W{u, vt = ^ [M\u, v) + L\u, v)] (18c) 



this becomes 



•^0 



a^iP) = — iu\ dpp/o(vp)e2 



:^^"P? 



(17b) 



As expected, these equations describe a circular dif- 
fraction pattern which is fully determined by a radial 
variable, c or v . The pattern is centered at the geometri- 
cal source image C, defined by c = v = 0, and a^^iP) is 
constant on any circle about C . The radius of the geo- 
metrically illuminated spot at the distance r from the 
aperture is a (ro + r)/ro, so that in the notation of Eq. 
(17a) the geometrical shadow boundary is defined by 
v = u. The parameter w, which relates the aperture ra- 
dius a to the wavelength A and the distances ro and r,^ 
can assume widely different values. For example, in the 
case of a classroom demonstration of Fresnel diffrac- 
tion, the parameters A = 500 nm, a = 0.\ mm, ro= r = 
100 mm are typical and in this case one has u = 0.87T. 
On the other hand, for limiting apertures used in a 
radiometer, parameters such as A = 500 nm, a = 5 mm, 
ro = r = 1 m are typical, and then one has u = 200tt. As 
will be shown later, the diffraction patterns encountered 
in these different cases are very different. For u ^ the 
diffraction pattern approaches the Fraunhofer limit 
(Airy function), and for w ^ co it approaches the limit 
of geometrical optics (rectangle function). (See Sec. 
4.2, Figs. 6a-d). 



^ It should be noted that, according to Eq. (13c), the distance r = z/ 
cosO depends on the location of the point of observation P so that, 
strictly speaking, u is not a constant if the diffraction pattern is ob- 
served at at fixed distance z from the aperture. This dependence of u 
on P is considered negligible in the Fresnel approximation. 



is the relative irradiance of the diffracted light and 

^^ X Im(Q;L) L(w, v) ..r... 

(Pi^(u , v) = arctan ^ ) , = - arctan ^ ", ' \ (1 8d) 
Re(aL) M(m,v) 



is the phase difference relative to the geometric field. 
The functions L(w, v) and M(w, v) appearing in these 
equations are defined by 



-L(m, v) = sm — + Vo(w, v)sm- - Vi(w, v)cos- 



= Ui(m, v)cos--i-U2(w, v)sin- , (19a) 



-M(w, v) = COS — + Vo(m, v)cos- - Vi(w, v)sm- 



= Ui(w, v)sin- - U2(w, v)cos-, (19b) 



where 



Vo(^/, v) = Jo(v) - Q'J2(v) + QVv) + - ..., 



(20a) 



V.(«,v)=Qj,(v)-Q3J3(v)-hQ%(v 



)+ 



(20b) 



Elsewhere in the literature, L(m, v) and M(m, v) are denoted by 
C(m, v) and S(m, v). This practice is not followed here in order to 
avoid confusion with the Fresnel integrals C{s) and S{s) of Eq. (8). 
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UK^,v) = r JiW- -h3(v)+ -%(v) + 



(20c) 



U2(u, v) = ( - )'J2(v) - I - ) Wv) + ( - )'J6(v) + 



(20d) 



Accordingly, 



Rc[ai^(u,u)]=-[l - Jo(m)]cos-, 



Im[Q;L(w, u)] = - 2^^ + Jo(w)]sin-, 



(23e) 



are Lommel functions of two variables, J„(v) being a 
Bessel function of the first kind and order n. 

For checking the accuracy of numerical results, it is 
useful to note the values of these expressions in special 

cases: 



|q;l(w,w)P = -[1 - 2Jo(w)cosw + J§(m)]. (23f) 



^(m , m ) = — arctan 



1 + Jo(m ) u 
1 — Jo(w) 2. 



(23g) 



a. In the limit m ^ 0, Lommel's equations simplify to 
the familiar Airy formula for Fraunhofer diffraction at a 
circular aperture. In this case, the entire diffraction pat- 
tern lies in the geometrical shadow and Eqs. (20c,d) are 
reduced to 

Ui(m,V) ^ Ji(v) U2(m,V) yQ ^ (p\ y ~ ^^hiv) 

U V ' U ' V ' 

(21) 

so that Airy's formula, U(P) ^ Ji(v)/v, is obtained from 
Eqs. (17a) and (3a,b). 

b. For V = 0, we have 



Re[Q;L(w , 0)] = ( 1 - cos - ) , Im[Q;L(w , 0)] = - sin - , 

(22a) 



Iq;l(w,0)P = 2 1 - cos- 



(p^(u , 0) = — arctan 



2) 
sm{u/2) 



_1 — cos(w/2) 
c. For V = w, the well-known relations [9] 



-cosw = - Jo(w) - h{u) + J4(w) ..., 



- sinw = Ji(m) — h(u) + h(u) ..., 



may be used to show that 



(22b) 
(22c) 



(23a) 
(23b) 



Vo(w, u) = - [Jo(w) + cosw], Vi(m, u) = - sinw, (23c) 



Ui(m, u) = -sinw, U2(w, u) = - [Jo(w) — cosm]. (23d) 



The use of Lommel's equations for numerical compu- 
tations is straightforward, provided that accurate values 
of the Bessel functions Jn(v) required for Eqs. (20a-d) 
are available and the convergence behavior of these 
equations is taken into consideration. 

When v/u or u/v are small, these expansions will 
converge on account of the monotonic decrease of (v/ 
m)" or (w/v)", provided of course that L and M are 
evaluated in terms of the Lommel functions Vo and Vi 
when V <u and in terms of Ui and U2 when v> u. 

When v/u or u/v are close to unity, Eqs. (20a-d) will 
converge on account of the relation J„ (v ) ^ when 
/2 ^^ 00. The manner in which this limit is approached is 
illustrated in Fig. 5. As n increases, J„(v) exhibits an 
oscillatory behavior before vanishing after passing 
through a pronounced final maximum at or below n = v. 
In this case, L and M can be evaluated in terms Vq, V\ 
or Ui,U2, although for consistency it is better to use the 
former method when v<u and the latter method when 
v> u. 

It follows that, for computations of the Lommel func- 
tions to m decimals, the expansions Eqs. (20a-d) can be 
truncated when either of the conditions, 




^■^55 



1P m 30 rl 

Fig. 5. Dependence of Bessel functions J„(v) on n. 



505 



Volume 103, Number 5, September-October 1998 

Journal of Research of the National Institute of Standards and Technology 



ij<^'<>- 



n > V and J„(v)<;r 10' 



(24a) 
(24b) 



are satisfied. 

The numerical results presented in this paper were 
obtained on a personal computer, using standard spread- 
sheet software^ (a 133 Mhz Pentium computer and Mi- 
crosoft Excel 7.0). It was found that this software pro- 
vides accurate values of the Bessel functions J„(v) 
needed for Eqs. (20a-d) without problems, but that the 
large number of them required to satisfy Eq. (24b) im- 
peded the speed of program execution when u is large 
and V ~ w . In addition, the capabilities of the personal 



computer were overtaxed by the fact that the diffraction 
patterns for large values of u are highly structured (see 
Figs. 6a-d), so that a very large number of data points 
had to be computed. For these reasons, Eqs. (18a) to 
(20d) were used in this work only for u ^ 300 while for 
larger values of the approximation of Sec. 4.3, below, 
was used. It should be emphasized that this limitation is 
unnecessary for larger computers. When sufficient com- 
puting power is available, the Lommel functions for 
large values of u can be evaluated efficiently by iterative 
use of recurrence relations for Bessel functions, begin- 
ning at the required large orders and iterating towards 
Jo(v) from above, as mentioned by Shirley and Datla 
[13]. Under these conditions, the fine structure of the 
diffraction pattern poses no difficulties. 




■It 



y^ 




2t 




1 v/u 1 .25 






tiS.i 



MW^^ 



m^<- 



-\- 



023 



5 75 



1 1^ 1.2S 



Fig. 6. Relative irradiances for circular apertures: a, b, c) Iq;l(m, v)I^ vs v/u for m = 1, 10, and 100 according to Sec. 4.2. d) IctschwC", v)l^ vs v/u 
u for u= 1000 according to Sec. 4.3. 



Certain commercial equipment, instruments, or materials are identi- 
fied in this paper to foster understanding. Such identification does not 
imply recommendation or endorsement by the National Institute of 
Standards and Technology, nor does it imply that the materials or 
equipment identified are necessarily the best available for the purpose. 
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The algorithm used in this work for u ^ 300 was as 
follows. 

1 . Define the value of m , a maximum value Vmax, a step 
size Av, and the desired decimal accuracy lO""". 

2. Let V = 0. Compute q;l(w, 0) from Eq. (22a). 

3. Let V = V + Av. Compute (u/2) Vo(m, v), (w/2)- 
Vi(m,v) from Eq. (20a,b), terminating when Eq. 
(24a or b) are satisfied. Compute L(m, v), M(m, v) 
from Eq. (19a,b) and q;l(m, v) from Eq. (18b). 

4. If v<M, go to Step 3. 

5. Compute q;l(w, u) from Eq. (23 e). 

6. Let V = V + Av. Compute (u/2) l]\(u, v), (u/2) 
U2(m,v) from Eq. (20c,d), terminating when Eq. 
(24a or b) are satisfied. Compute h(u, v), M(m, v) 
from Eq. (19a,b) and ai^(u, v) from Eq. (18b). 

7. If V < Vmax, go to Step 6. Else, stop. 

The relative irradiances Iq;l(w, v)I^ computed in this 
manner for w = 1, 10, and 100 are plotted in Figs. 6a-c. 

4.3 Schwarzschilds's Approximation 

The above-mentioned computational problems en- 
countered with Lommel's solution when u is large can 
be avoided by using an asymptotic approximation for 
ap(P) derived by Schwarzschild [8] in a paper on dif- 
fraction effects in defocused telescopes. As this paper is 
no longer readily available, its contents will be outlined 
here. 

Schwarzschild considered the integral 



W 



lu 

''2^ 



= £JoHo 

f dp... - J 



d;^e 



■i(ufP'/2-vpcosx+y^^'^u) 



= W,- W2. 



(25a) 



so that his approximation of Eq. (17b) is given by'* 



ir , 



q;l(m, v) ~ aschw(u, v) = eiu W^ 



(25b) 



The integral Wi is readily shown to be equal to 

W, = \, (26a) 

and by evaluating the ;^-integral of W2 as in Eq. (16b) 
and then substituting the asymptotic expression [9] 



^° It is well known that the wave functions for diffraction with and 
without a lens are complex conjugates. 



Jo(vp) = y— cos(vp - Jj, vp » 1 



(26b) 



one obtains 



V277V L J\ 



-{p^-vluY 



+ e 4 iM 



Letting 



/= 



p + v/u ' 



g=- C 2 



v/uY 



{f^vluY 



(26c) 



(27a) 



the first integral in Eq. (26c) can be expressed in the 
form // d^ so that 



i(«+v)^ 
111 



.^--(^v/«)2_ 



■m\ dp\/^e->^^«)' = ^ 
Ji 1+v 

"^J, \p + v/u) 
The second integral in Eq. (26b) can be written as 

'm\ dpV^e-f^^-^'"^' 
=m\ dp(Vp - Vv7^)e-f^^-^'"^' 



(27b) 



iVwj 



dpe 



-f(p-v/«)2 



(28a) 



where the first term can again be evaluated by partial 
integration, using 



/= 



1 



Vp — vv/u _ 
p-v/u " V^ + Vw^ 



(28b) 



and the second term is a complex Fresnel integral [Eq. 
(8)]. In this manner, Schwarzschild found 

iu\ dpV^e-^t^-^^">' 

= iV77v{F*(oo) - F*[V^(1 - v/u)]} 

+ ^^ + I df ^-^ )e-f*''->l (28c) 
1 + Vv/u J 1 \ Vp + Vv/u/ 

He noted that, by further partial integrations, the result- 
ing expression for W2 would become an asymptotic ex- 
pansion in negative powers of u but that there was no 
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point in taking the trouble as the last terms in Eqs. (27b) 
and (28c) are already negligibly small for practical pur- 
poses. Therefore, 



1 



W= 1 -W2-^:r + ^e-TF*[VM/7r(l -vlu)] 

2 V2 



[jU+Vr ITT i(« — V) JTT 

e 2h "'"4 e 2h 4 



1 + vlu 1 + - 



u,v» 1. (29) 



Schwarzschild estimated that this expression is accurate 
to 0.005 if M = 100, vlu > 0.2; or w = 300, vlu > 30. 

When put in the form of Eq. (25b), Schwarzschild's 
approximation becomes 



«f(m , v) ~ Q;schw(M , v) = - e ^* - 



v^' 



-i(S-T7/4) 



F(^) 



1+v/m TTVv/mX^' 



v»l, (30a) 



Re[a;schw(w, v)]= - cosS ^ [C(^)sin(5 - ttM) 

2 V2 



5(^)C0S(5 - 7T/4)] 



COS(j6- + 77/4) 



COS(j6+ - 77/4) 

1 + v/m 



, M , V » 1 , 



(30b) 



Im[Q;schw(M, v)]= - - sinS j= [C(^)cos(6 - 77/4) 

2 V2 



+ 5(^)sin(5-77/4)] -■ 



sin(j6.H - 77/4) 
1 + v/m 



sin(j6- + 77/4) 



, M, V» 1.., 



(30c) 



where 



5 = |-,^^=^±v,^ = V^(l - vlu). (30d) 



The use of these expressions on a computer is simple. 
The only caveat is that they are not valid for small values 
of vlu so that Lommel's equations must still be used 
below a suitably chosen minimum value v = Vmin- The 
choice of Vmin can be based on the following table, which 
was obtained by computing the residuals 
A = laschwl^ — \oiif for selected values of u. 



For such small values of vlu , only a few terms of Eqs. 
(20a,b) are sufficient to obtain ckl with comparable ac- 
curacy. Thus, the following procedure will provide the 
entire diffraction pattern. 

1. Define m, v^^^, v^^x, and Av. 

2. Let V = 0. Compute q;l(0, v) from Eq, (22a). 

3. Let V = V + Av. Compute (m/2)-Vo(m, v), (m/2)- 
Vi(m,v) from Eqs. (20a,b), terminating after the 
third terms. Compute L(m,v), M(m,v) from Eqs. 
(19a,b) and q;l(w, v) from Eq. (18b). 

4. If V < Vmin, go to Step 3. 

5. Let V = V + Av. Compute 8, (5±, s from Eq. (30d). 
Use the algorithm of Sec. 3.1 to find C{s), S(s). 
Compute aschw(u, v) from Eqs. (30b,c). 

6. If V < Vmax» go to Step 5. Else, stop. 

The values of lQ;schw(w, v)P computed by this algorithm 
for u = 1000 are shown in Fig. 6d. They were found to 
be in excellent agreement with the values of Iq;l(m, v)P 
obtained from Lommel's solution. ^^ 

It should be emphasized that Schwarzschild's approx- 
imation is different from a superficially similar but sig- 
nificantly less accurate asymptotic approximation of the 
diffraction integral (17b) cited by Focke [7] and used by 
Blevin [14], Steel, De, and Bell [16], and Boivin [16] in 
their work on diffraction errors in radiometry. The re- 
spective accuracies of the Schwarzschild and Focke ap- 
proximations can be assessed by comparing the relative 
irradiances at the shadow boundary,^^ 

\c^Schw(u, m)P = - 1 — ^J COSMCOsI M — — j 



2 

+ COS 

77M 



("-f). 



(31a) 



lQ;Focke(M, m)I^ = 7 1 7= (cos2m + sin2M) -\-- — , 

4 L V77M 277MJ 

(31b) 

to the exact value Iq;l(w, m)P given by Eq. (23f). From 
the ratios plotted in Fig. 7 it follows that the 
Schwarzschild values are accurate to 0.1 % for w > 10, 
while for w = 100 the Focke values are still off by 6 %. 



^^ This comparison was kindly performed by Dr. Eric Shirley of NIST. 
^^ Equations (3 1 a,b) follow from Eq. (30a) of this paper and Eq. (7) of 
Ref. [13], respectively. 



u 


A < 0.01 


A < 0.001 


30 


if v/u ^ 0.23 


if v/u ^ 0.7 


100 


if v/u ^ 0.06 


if v/u ^ 0.27 
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Fig. 7. Irradiance ratios, lQ:schw(w, w)l /\ai^(u^ u)\ and lQ:Fockc(«* w)l / 
Iq;l(m, m)I , vs u. 
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